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ABSTRACT 

We construct models of triaxial galactic nuclei containing central black holes using the method of orbital super- 
position, then verify their stability by advancing A^-body realizations of the models forward in time. We assume 
a power-law form for the stellar density, p oc r~ 7 , with 7=1 and 7 = 2; these values correspond approximately to 
the nuclear density profiles of bright and faint galaxies respectively. Equidensity surfaces are ellipsoids with fixed 
axis ratios. The central black hole is represented by a Newtonian point mass. We consider three triaxial shapes for 
each value of 7: almost prolate, almost oblate and maximally triaxial. Two kinds of orbital solution are attempted 
for each mass model: the first including only regular orbits, the second including chaotic orbits as well. We find 
that stable configurations exist, for both values of 7, in the maximally triaxial and nearly-oblate cases; however 
steady-state solutions in the nearly-prolate geometry could not be found. A large fraction of the mass, of order 
50% or more, could be assigned to the chaotic orbits without inducing evolution. Our results demonstrate that 
triaxiality may persist even within the sphere of influence of the central black hole, and that chaotic orbits may 
constitute an important building block of galactic nuclei. 

Keywords: galaxies: elliptical and lenticular — galaxies: structure — galaxies: nuclei — stellar dynamics 



1. INTRODUCTION 

Schwarzschild (1979) demonstrated how to construct self- 
consistent models of stellar systems in the absence of analytic 
expressions for the orbital integrals. His method consists of 
three steps, i) Represent the stellar system by a smooth den- 
sity law and divide it into discrete cells; ii) compute a library 
of orbits in the potential corresponding to the assumed density 
law, and record the time spent by each orbit in the cells; iii) find 
a linear combination of orbits that reproduces the cell masses. 
Using his method, Schwarzschild (1979, 1982) demonstrated 
self-consistency of triaxial mass models with and without fig- 
ure rotation. Most of the orbits in his solutions were regular, i.e. 
non-chaotic (Merritt 1980). Subsequently, Statler (1987) found 
a variety of self-consistent solutions for the integrable, or "per- 
fect," triaxial mass models in which all orbits are regular. 

Models like these with large, constant-density cores are now 
known to be poor representations of ellipti cal galaxies, almost 
all of which have high central densities ( prane et al. 1993 ; 
Ferrarese et al. 1994). Stellar densities rise toward the cen- 
ter approximately as power laws, p oc r~ 7 . Fainter galaxies 
have steeper cusps, 7 w 2, while brighter galaxies have weaker 
cusps, < 7 < 1, and exhibit an obvious break in the surface 
brightness profile. Following this discovery, Schwarzschild 
(1993) investigated triaxial models with singular density pro- 
files, p ~ r" 2 , and Merritt & Fridman (1996) constructed self- 
consistent solutions for triaxial galaxies with both weak and 
strong central cusps. A significant portion of the phase space 
in these models was found to be occupied by stochastic orbits. 
Furthermore, triaxial self-consistency could sometimes only be 
achieved by including some stochastic orbits. Models contain- 
ing stochastic orbits can represent bona-fide equilibria as long 
as the stochastic orbits are represented as fully-mixed ensem- 
bles (Merritt & Fridman 1996; Merritt & Valluri 1996). 

In the last decade, evidence has grown that supermassive 
black holes are generic components of galactic nuclei (Ho 
1999). There are roughly a dozen galaxies in which a com- 



pelling case for the presence of a supermassive black hole can 
be made based on the kinematics of stars or gas (Merritt & Fer- 
rarese 2001), as well as a number of active galactic nuclei in 
which the kinematics of the broad emission line region imply 
the existence of a supermassive black hole (Peterson 2002). In- 
ferred masses range from ~ 10 6 Alo to ~ \Q 95 M.q and corre- 
late well with stellar velocity dispersions (e.g. Ferrarese et al. 



2001) and bulge luminosities (e.g. |McLure & Dunlop 2002 ). 

The possibility of maintaining triaxiality within a galactic nu- 
cleus containing a supermassive black hole remains a topic of 
interest. Very close to the black hole, the gravitational force 
can be considered a perturbation to t he Kepler problem and 
the phase space is essenti a lly regular (|Merritt & Valluri 1999 



Sambhus & Sridhar 2000j ; [Poon & Merritt 200 1| , hereafter Pa- 
per I). Farther from the black hole, the fraction of chaotic orbits 
increases, up to a radius where the enclosed stellar mass is a 
few times the black hole mass; beyond this radius essentially 
all centrophilic orbits are chaotic (Paper I). The tube orbits re- 
main mostly regular since they avoid the destabilizing center. 
The persistence of regular orbits throughout the region where 
the gravitational force from the black hole dominates, leaves 
open the possibility of constructing self-consistent solutions. 
Furthermore there is growing observational evidence for the ex- 
ist ence of bar-like distort ions at the very centers of galaxies (e. 
g. |Erwin & Sparke 2002|). 



In Paper II (Poon & Merritt 2002), we presented preliminary 
results showing that self-consistent and stable triaxial equilib- 
ria could be constructed for power-law nuclei with certain axis 
ratios. In this paper, we present a more detailed investigation of 
triaxial black-hole nuclei. We find that stationary solutions are 
possible only for certain shapes; mass models that are too near 
to prolate axisymmetry always evolve toward axisymmetry. 

The properties of the mass models are presented in §2. Or- 
bital solutions for various shapes and density profiles are pre- 
sented in §3 and their stability is tested by N-body simulation 
in §4 and §5. Limits on the chaotic mass fraction are discussed 



2 



in §6. §7 discusses some implications for the nuclear dynamics 
of galaxies. 



that 



2. MASS MODEL 

We model the stellar distribution by the density law: 



= p a m 



b 2 



(1) 

(2) 



The equidensity surfaces are concentric ellipsoids with fixed 
axis ratios a : b : c and the radial profile is a power law with 
index -7. We define the outer surface by the ellipsoid m = m ou , 
and measure the triaxiality via the index T, where: 



T = 



(3) 



Oblate and prolate galaxies have T = and T = 1 respectively. 
T = 0.5 corresponds to a "maximally triaxial" nucleus. 

We consider two types of nuclei, the weak cusp (7=1) and 
the strong cusp (7 = 2), which correspond roughly to the den- 
sity profiles observed at the centers of bright and faint elliptical 
galaxies respectively. For each value of 7, we consider three 
shapes: almost oblate (T = 0.25, a : b : c= 1.0 : 0.9 : 0.5); max- 
imally triaxial (T = 0.50, a : b : c= 1.0 : 0.79 : 0.5); and almost 
prolate (T = 0.75, a : b : c = 1.0 : 0.66 : 0.5). All models have 
c/a = 0.5. 

The black hole is represented by a central point mass with 
Mhh = 1, which imposes a scale to the otherwise scale-free stel- 
lar mass model. We define two characteristic radii associated 
with the presence of the black hole (Tables 1, 2). i) r g is de- 
fined such that the enclosed stellar mass within an ellipsoid with 
m = r g is equal to that of the black hole. For T = 0.5, r g = 0.64 
for 7=1 and r g = 0.20 for 7 = 2. ii) r c /, is the radius beyond 
which the regular, boxlike orbits become almost all stochastic. 
For T = 0.50, r c h ~ 2r g for 7=1 and r c h ~ 6r g for 7 = 2. 

The gravitational potential can be ob tained from Chan- 
drasekhar's theorem ( JChandrasekhar 1969| ), which is valid for 
density laws that are stratified on similar ellipsoids. The corre- 
sponding forces may be obtained in analytical form by taking 
partial derivatives of the potential. The details are given in §2 
of Paper I. 

Our orbital solutions require a finite outer radius. Our aim 
was to explore the possibility of maintaining triaxiality out to a 
radius of at least r c /,; hence we chose the outer surface of our 
model to be large enough that almost all of the density at r C h in 
a real galaxy would be contributed by orbits with apocenters r+ 
lying below this surface. 

To estimate m ollt , we considered a spherical galaxy with den- 
sity 



(4) 



The isotropic distribution function corresponding to the density 
law (0) can be found using Eddington's formula, 
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with $ the potential generated by p*(r) and by the central point 
mass, and = linv^oo $(r). We then change variables such 



dp 



sph 
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r 2 v r 



d(E,L 2 ) 



0(r+,r_) 



dr + dr_, 



(6) 



with v r = ^/2(E-L 2 /2r 2 - <&(r) the radial velocity, and r+ and 
r_ the apocenter and pericenter of an orbit with energy E and 
angular momentum L. The Jacobian is 



d(E,L 2 ) 



According to (|j), we define 

2tt/; 
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u(r+,r) : 



d(E,L 2 ) 



dr. 



p7 h (r) Jo 



d{r+,r.) 
g(r+,r)dr+. 



(8) 
(9) 



u(r+, r) is the fraction of the mass at r contributed by orbits with 
apocenters r < r + . 



7=1 



u(r + , r) 



u(r + , r) 




7=2 




12 3 4 

r + /r 

FIG. 1 . — u(r+, r) as a function of r+/r for 7=1 and 2 and for different radii, 
as indicated. 



Figure |lJplots u(r+, r) as a function of r+/ r for both weak and 
strong cusp nuclei. The red curves correspond to infinite r, i.e. 
no influence from the black hole. The black curves correspond 
to very small r, where the potential is almost Keplerian. u(r +1 r) 
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for the weak cusp case is smaller than that for the strong cusp 
at a given r + /r, because the density of the strong cusp falls off 
faster, thus the mass at a given r has to rely more heavily on 
nearby orbits. 

Based on these results, we took m out to be 5r c h for 7=1 and 
3r c h for 7 = 2, so that about 70%(75%) of the mass at r c h is 
accounted for in the weak (strong) cusp cases respectively. Ta- 
bles |l]and I give the values of the characteristic radii in model 
units. 

Table 1 
Model Parameters for 7 = 1 





T = 0.25 


r = o.50 


T = 0.75 




0.594 


0.635 


0.694 


r c h 


1.057 


1.128 


1.234 


m out 


5.284 


5.642 


6.168 






Table 2 






Model Parameters for 7 = 2 




T = 0.25 


T = 0.50 


7 = 0.75 


r g 


0.177 


0.201 


0.241 


r c h 


1.114 


1.270 


1.518 


m out 


3.342 


3.811 


4.555 



3. CONSTRUCTION OF ORBITAL SOLUTIONS 

We followed standard procedures for constructing the 
Schwarzschild solutions. The cells used to define the mass dis- 
tribution were defined as follows. Each model with outer sur- 
face m ou , was divided into 64 shells. The inner 63 shells were 
equidensity surfaces; the outermost shell was an equipotential 
surface, in order to accomodate chaotic orbits which fill regions 
defined by equipotential surfaces. Shells were more closely 
spaced near the center. The 42nd shell in each model corre- 
sponded to r c h- Each shell was further divided into 48 angular 
cells per octant as Merritt & Fridman (1996), giving a total of 
3072 cells per octant. Due to the symmetry of the problem, it 
is necessary to consider only a single octant when constructing 
the orbital solutions. 

Orbits were computed in two initial condition spaces: sta- 
tionary start space, which yields mostly centrophilic orbits 
(pyramids, stochastic orbits); and X - Z start space, which 
yields mostly tube orbits. Orbital energies were selected from 
a grid of 42 (52) values for 7=1 (2), defined as the energies 
of equipotential surfaces which were spaced similarly in radius 
to the equidensity shells. The outermost energy shell, which is 
also an equipotential surface, intersects the .x-axis at x = m out 
for both mass models. X—Z start space consists of orbits which 
begin on the X-Z plane with v x = v, = 0; stationary start space 
consists of orbits which begin with zero velocity. A total of 
18144 (22464) orbits were integrated for 100 dynamical times 
for 7 = 1(2) and their contributions to the masses in the cells 
were recorded. In order to distinguish regular from stochastic 
trajectories, the largest Liapounov exponent was computed for 
each orbit. 

Twelve equations - six representing the unperturbed motion 
and six the linearized perturbations - were integrated using the 
routine RADAU of Hairer (1996), which is a variable time step, 
implicit Runge-Kutta Scheme which automatically switches be- 
tween orders of 5, 9, and 13. Energy was typically conserved 
to a few parts in 10 9 over 100 orbital periods. 



An orbit was considered regular if the largest Liapunov expo- 
nent A satisfied A7b < 10~° 9 . The dynamical time 7b is defined 
as the period of a circular orbit of the same energy in the equiv- 
alent spherical potential, which is defined to have a scale length 
(abcp (Paper I). This threshold was determined empirically by 
making histograms of Liapounov exponents of mono-energetic 
orbits at various energies; A7b as 10 -0 9 was found to always 
separate the two peaks in the histogram corresponding to regu- 
lar and chaotic orbits. A large fraction of the computed orbits 
were found to be stochastic, as shown in Table 0. 

Table 3 

Fraction of regular orbits in the orbit libraries 





7=1 


7 = 2 


T = 0.25 


0.58 


0.72 


T = 0.50 


0.54 


0.67 


T = 0.75 


0.49 


0.65 



We then found the linear combination of orbits which best 
reproduced the cell masses. We did this by varying the orbital 
occupation numbers C, to minimize the quantity \ 2 , defined as: 

N M 

X 2 = ^(£>,-$> ; ,C,) 2 . (10) 

l=\ i=\ 

Bn is the time spent by the j-th orbit in the Z-th cell, D\ is the 
mass of the Z-th cell, and C-, is the occupation number of the 
Z-th orbit. The quadratic programming problem was solved us- 
ing the NAG Fortran library routine E04NFC. We measured the 
discrepancies of the orbital solutions by the parameter: 

i=\ 1 >=i 

the mean error in the cell masses. 

We constructed two solutions for each mass model: one 
using only regular orbits, and the other using both regular 
and chaotic orbits. Schwarzschild (1993) was the first to in- 
clude chaotic orbits in self-consistent solutions. Schwarzschild 
treated the chaotic orbits like regular orbits, giving each chaotic 
orbit its own occupation number, even though many of the 
chaotic orbits in his models were "sticky" and did not reach 
a time-averaged steady state during the interval of integration. 
He justified this practice by re-integrating all chaotic orbits with 
non-zero occupation numbers for longer intervals in the fixed 
potential; the extreme error in the cell masses was found to in- 
crease from the original ~ 1% to ~ 10%. Merritt & Fridman 
(1996) noted that fully-mixed chaotic orbits, i.e. chaotic orbits 
that uniformly fill their accessible phase space region, are bona- 
fide building blocks for steady-state galaxies and approximated 
such building blocks by constructing ensemble averages of the 
chaotic orbits at each energy. 

We followed Schwarzschild in allowing different chaotic or- 
bits at a given energy to have different occupation numbers. In 
our potentials, the chaotic orbits were observed to very quickly 
fill the region accessible to them: their "mixing times" were 
generally much shorter than the integration interval. Hence we 
expected that our chaotic orbits would individually constitute 
bona fide building blocks for a stationary solution, without the 
need to construct ensemble averages as in Merritt & Fridman 
(1996). This expectation was confirmed by the ZV-body integra- 
tions described below. 
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FIG. 2. — Schematic diagram showing the three kinds of shells used for 
constraining the orbital solutions, (i) The innermost shells with all angular 
constraints imposed, (ii) Intermediate shells with angular constraints ignored, 
(iii) The outermost shells, which are ignored. The highest-energy shell is an 
equipotential surface while the other shells are equidensity ellipsoids. 

We could not find quadratic -programming solutions which 
exactly reproduced the cell masses in all of the cells, includ- 
ing the outermost ones. This is probably because the number of 
orbits visiting the outer shells is small, giving the quadratic pro- 
gramming algorithm little freedom to fit the cell masses there. 
Indeed in all of our self-consistent solutions, most of the con- 
tribution to A 2 came from the outermost cells. We therefore 
considered solutions in which the outermost cells were treated 
in various ways. Figure [l] illustrates the idea. We defined three 
kinds of mass shells: (i) the innermost shells with all the con- 
straints imposed, i.e. all of the angular cells included; (ii) in- 
termediate shells for which only the total mass was fit, with 
angular details ignored; (iii) the outermost shells which were 
excluded from the fit. We carried out numerous tests in which 
we attempted to fit various combinations of constraints. As ex- 
pected, the discrepancy A defined by the innermost shells al- 
ways decreased as the total number of constraints decreased. 
For instance, for the weak cusp model with T = 0.50, including 
all the orbits allowed the innermost 60 shells (r ~ 4r c /,) to be fit 
to machine precision when the outermost four shells were ig- 
nored. When these "exact" solutions were advanced forward in 
time, however, they were found to exhibit significant evolution 
at large radii, due presumably to the poor fit in the outermost 
cells. We discuss this further in §5. Experiments like this per- 
suaded us to focus on solutions that were not "exact" but rather 
were constrained to reproduce the densities in all cells within 
m out to as high an accuracy as possible. Such solutions exhib- 
ited smaller fractional errors in the inner cells (r < r c h) than 
in the outer ones. Models constructed in this way were found 
to evolve much less than the "exact" solutions and provide the 
basis for the discussion below. 

The orbital content and the precision A of the solutions are 
represented in various ways in Tables and in Figures ||- 
||. We classify orbits into one of four families (cf. Paper I). 
x- and z-tubes are regular orbits that circulate around the long 
and short axes of the figure, respectively. Pyramid orbits are 



the closest analogs to box orbits in these potentials; they can 
be described as eccentric Keplerian ellipses that precess due to 
torques from the stellar potential. Their major elongation is 
contrary to that of the figure. We will indiscriminately refer to 
these orbits as "pyramids" and "box orbits" in what follows. 
The fourth family consists of the chaotic orbits. 

Tables [i] - 1| give the mass fractions within various radii con- 
tributed by the four types of orbit, z-tube orbits (short-axis 
tubes) are the biggest contributors in all of the regular-orbit so- 
lutions, making up at least 50% of the total mass and as much as 
80% in the nearly-oblate (T = 0.25) models. Pyramid orbits are 
next in importance; their contribution reaches ~ 40% for 7=1 
and ~ 15% for 7 = 2. jc-tubes (long-axis tubes) are relatively 
unimportant in the models with T = 0.25 and 0.5, making up 
only a few percent of the total mass. In the nearly prolate mod- 
els (T = 0.75), their contribution increases to ~ 20%, however 
we argue below that these prolate solutions do not represent true 
equilibria. 

When chaotic orbits are included in the orbit libraries, the 
character of the solutions changes substantially: at least 40% 
of the mass is assigned to chaotic orbits by the quadratic pro- 
gramming algorithm, and as much as 60% in the models with 
7=1 and T = (0.25,0.5). The inclusion of chaotic orbits low- 
ers the mean error A by a modest amount in both the weak and 
strong cusp cases (Tables 4-6). Much of the mass assigned to 
the chaotic orbits appears to be "taken" from the pyramid or- 
bits, as expected. To our knowledge, these solutions contain a 
larger fraction of chaotic orbits than in any other self-consistent 
galaxy models. Since the time-averaged shapes of the chaotic 
orbits are nearly spherical, the existence of these solutions im- 
plies that the regular orbits have more than enough freedom to 
reproduce the triaxial shapes, and can do so even when forced 
to compensate for the inauspicious shapes of the chaotic orbits. 

Another way to represent the orbital makeup of the solutions 
is via the distribution of orbital energies. Let Mj(E)dE be the 
mass in orbits from the ith orbital family whose energies lie 
in the range E to E + dE. The cumulative mass is given by 

M,-(< E) = J^ oo M i (E)dE. Figure | shows M,(< E) for the self- 
consistent solutions with T = 0.5. We show for comparison 
M sph (< E) computed for the equivalent spherical models de- 
fined in the Appendix. One expects that M sph (< E) w J2i M '(< 
E), and Figure [3] verifies that this is correct. At high energies, 
there are discrepancies between the Schwarzschild solutions 
and the predictions of the spherical models. This is because 
the Schwarzschild solutions are truncated; thus, orbits with the 
highest energies are forced to make large contributions, which 
in turn lowers the contributions from orbits at lower energies. 
This effect is more significant for the weak cusp models which 
place most of their mass at large radii. Nevertheless, at low en- 
ergies even the 7=1 solutions have nearly the same M(< E) 
dependence as the equivalent spherical models. 
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o 10 20 -30 -20 -10 

E E 

FIG. 3. — Cumulative energy distributions of the various orbital families, in self-consistent solutions with T = 0.5. The symbols "B", "X", "Z", "C" denote the 
mass contributed by box, ;t-tube, z-tube and chaotic orbits respectively; black dots give the total. Solid (dotted) lines show M(< E) for the equivalent spherical 
models defined in the Appendix, with (without) the central black hole. 



Table 4 



Orbital Content of Solutions with T = 0.25 





z-tubes 


x-tubes 


pyramids 


chaotic 


7=1: (regular) 
log(A) = -1.115 










r <0.5 


0.70 


0.06 


0.24 




r < 1.0 


0.62 


0.06 


0.33 




r < 1.5 


0.59 


0.05 


0.36 




7=1: (all) 
log(A) = -1.318 
r <0.5 


0.40 


0.03 


0.07 


0.50 


r < 1.0 


0.34 


0.02 


0.06 


0.59 


r < 1.5 


0.30 


0.01 


0.06 


0.63 


7 = 2: (regular) 
log(A) = -1.437 
r <0.4 


0.81 


0.07 


0.11 




r <0.8 


0.80 


0.08 


0.12 




r < 1.2 


0.80 


0.08 


0.12 




7 = 2: (all) 
log(A) = -2.633 
r <0.4 


0.56 


0.03 


0.07 


0.34 


r <0.8 


0.55 


0.03 


0.07 


0.35 


r < 1.2 


0.53 


0.03 


0.07 


0.37 
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Table 5 

Orbital Content of Solutions with T = 0.50 





z-tubes 


x-tubes 


pyramids 


chaotic 


7=1: (regular) 










log(A) = -1.169 










r < 0.5 


0.64 


0.13 


0.23 


— 


r < 1.0 


0.56 


0.08 


0.36 


— 


r < 1.5 


0.52 


0.07 


0.41 


— 


7=1: (all) 










log(A) = -1.307 










r < 0.5 


0.29 


0.07 


0.11 


0.54 


r < 1.0 


0.27 


0.03 


0.10 


0.60 


r < 1.5 


0.25 


0.02 


0.10 


0.62 


7 = 2: (regular) 










log(A) = -1.295 










r <0.4 


0.73 


0.11 


0.16 


— 


r <0.8 


0.73 


0.11 


0.15 




r < 1.2 


0.73 


0.12 


0.15 




7 = 2: (all) 










log(A) = -2.567 










r <0.4 


0.46 


0.05 


0.09 


0.40 


r <0.8 


0.45 


0.05 


0.06 


0.44 


r < 1.2 


0.44 


0.05 


0.05 


0.46 




Table 6 






Orbital Content of Solutions with T = 


0.75 




z-tubes 


x-tubes 


pyramids 


chaotic 


7=1: (regular) 










log(A) = -1.379 










r <0.5 


0.44 


0.27 


0.29 


— 


r< 1.0 


0.48 


0.20 


0.32 


— 


r < 1.5 


0.47 


0.16 


0.37 


— 


7=1: (all) 










log(A) = -1.248 










r <0.5 


0.27 


0.14 


0.16 


0.42 


r < 1.0 


0.26 


0.07 


0.14 


0.53 


r < 1.5 


0.22 


0.06 


0.16 


0.56 


7 = 2: (regular) 










log(A) = -1.221 










r <0.4 


0.57 


0.20 


0.23 




r <0.8 


0.59 


0.21 


0.20 




r < 1.2 


0.58 


0.22 


0.20 




7 = 2: (all) 










log(A) = -2.471 










r <0.4 


0.38 


0.13 


0.11 


0.38 


r <0.8 


0.35 


0.12 


0.09 


0.44 


r < 1.2 


0.33 


0.13 


0.08 


0.46 



Figure [I] presents yet another representation of the orbital 
populations of our self-consistent solutions. We plot the cumu- 
lative mass fraction F contributed by different kinds of orbits 
to different shells of the model. Box-orbits are blue, tube orbits 
are orange and chaotic orbits are purple; higher-energy orbits 
are represented by darker shades. Color bars relate the shades to 
the energy of the orbits. Numbers below the color bar indicate 
radii where the equidensity shells and equipotentials intersect 
the x-axis; note that the scale is not linear due to the higher reso- 
lution at lower radii. For instance, in the weak cusp model with 
T = 0.50, shell 20 is an equidensity ellipsoid with x-intercept 
0.78, and regular box orbits with energy E = $(0.78,0.0,0.0) 
are represented by the lightest blue color. 



At low energies, many of the z-tube orbits appear as saucers, 
i.e. 2 : 1 resonant orbits; at high energies, higher-order reso- 
nances are important for these orbits. It is surprising that x- 
tube orbits do not dominate the nearly-prolate solutions with 
T = 0.75. However configuration-space plots of the x-tube or- 
bits show that almost all of them are elongated contrary to the 
prolate figure: they are thin circular rings lying near the y— z 
plane. Even in the nearly-prolate solutions, most of the mass 
is contributed by z-tube orbits; the remaining contributions are 
mostly from high-energy box orbits, many of them associated 
with resonances. In the nearly-prolate solutions, the most im- 
portant resonances are the "fish" (2 : 3) and the "pretzels" (3:4). 
In the nearly-oblate solutions, the fish and the "banana" (2 : 1) 
resonances dominate. Many of the z-tubes and the high-energy 
boxes are replaced after the introduction of chaotic orbits (Fig- 
ure ||). 

As discussed in more detail below, we were not able to pre- 
dict the long-term stability of a model based on its value of 
A. Some models with fairly large A's were found to exhibit al- 
most no evolution, while other models with smaller A's evolved 
significantly. This suggests that the results of Schwarzschild 
modelling should be interpreted with caution in cases where 
the long-term stability of the solution has not been tested. 

4. W-BODY MODELS 

Each of the solutions described above was found via a min- 
imization of the quantity \ 2 describing the sum of the squared 
errors in the cell masses (equation 10). While the magnitude of 
X 2 might be expected to correlate with the quality of the solu- 
tion, there is really no way to know how small \ 2 must be in or- 
der for the solution to represent a bona fide steady state, or how 
large a value of \ 2 is consistent with the existence of a smooth 
equilibrium solution. One could require that each of the cell 
masses be fit exactly; but we were never able to do this when 
the outermost shells were included, and in any case, the dis- 
crete representation of the density renders the interpretation of 
an "exact" solution problematic. There are other uncertainties 
as well; for instance, both "sticky" chaotic orbits and nearly- 
resonant, regular orbits may require very long integration times 
before their cell masses approximate the steady-state values. 

We tested whether our orbital superpositions represent true 
equilibria by realizing them as A^-body models and integrating 
them forward in time in the gravitational potential computed 
from the N bodies themselves, and from a point mass represent- 
ing the black hole. We prepared initial conditions correspond- 
ing to each of the orbital solutions by re-integrating orbits with 
non-zero occupation numbers and storing their positions and 
velocities at fixed time intervals. The sense of rotation of the 
tube orbits was chosen randomly. We used ~2x 10 6 particles 
for representing the weak-cusp models and ~2x 10 5 particles 
for the strong-cusp models; a smaller N was chosen for 7 = 2 
since the A^-body integrations were slower for the more con- 
densed model. We int egrated the initial conditions forw ard in 
time using GADGET ( |Springel, Yoshida & White 200 1[ ), a par- 
allel tree code with variable time steps. The particle represent- 
ing the black hole was allowed to move in response to the forces 
from the "stars." The softening length was 0.005(0.003) for 
the weak (strong) cusp cases. Energy was conserved to within 
~ 0.5% for the weak cusp models and ~ 1 % for the strong cusp 
models. 

As our primary index of evolution, we computed the axis ra- 
tios of the models as a function of radius using the iterative 
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FIG. 4. — Cumulative mass fraction F contributed by different kinds of orbits to different shells of the triaxial solutions. Box orbits are blue, tube orbits (both z- 
and .v-tubes) are orange, and chaotic orbits are purple. Higher energies are represented by darker shades, defined according to the .v-intercept of the equipotential 
that has the same energy as the orbit. Numbers below the color bar indicate radii where the equidensity shells intersect the x-axis. Figures labelled "reg" respresent 
solutions constructed only using regular orbits. 



procedure described by Dubinski & Carlberg (1991), as fol- 
lows, (i) The moment of inertia matrix of the particles enclosed 
by a sphere of radius r is calculated, (ii) Axis ratios are as- 
signed as a = yjm u /m max , b = y / m 2 2/m mu , c = y/ 'm 33 / 'm max , 
where m,,- are the principal moments of inertia and m max = 



maxfmn, JM22, WZ33}. (iii) Particles are enclosed within the el- 
lipsoid x 2 /a 2 +y 2 /b + z 2 /c 2 = r 2 and step (ii) is iterated, until 
the axis ratios converge. 

Figures |5| and ^ show the evolution of the axis ratios and den- 
sity profiles of the various solutions. None of the density pro- 
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files showed signficant evolution: the radial distribution of mass 
did not evolve for any of the models. However some of the so- 
lutions showed significant evolution in their shapes. For T = 0.5 
(maximum triaxiality) and 7 = ( 1 , 2), we observed no significant 
evolution in the axis ratios, either for the regular-orbit solutions 



or the solutions containing chaotic orbits. We conclude that 
the maximally-triaxial solutions represent bona-fide equilibria. 
For the weak-cusp model with T = 0.25 (nearly oblate), there 
is some evolution at late times in the large-radius axis ratios, 
while the nearly-oblate model with 7 = 2 hardly evolves. Con- 
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FIG. 6. — Evolution of density profiles for the orbital solutions of Figure m Different colored curves correspond to the times shown in the upper left-hand panel, 
in units of the total integration time. Dotted lines are the density profile of trie underlying mass model. 



tour plots of the particle distribution suggest that the evolution 
for 7 = 1 is due to the relatively large errors in the cell masses 
at large radii; as time goes on these errors propagate to smaller 
radii. In spite of these fluctuations, however, we judged both of 
the T = 0.25 solutions to be stable. 

By contrast, all solutions with T = 0.75 (nearly prolate) were 
found to exhibit substantial evolution, reaching nearly axisym- 
metric (oblate) shapes by the final time step. (We note that 
for the weak-cusp solution with T = 0.75, even the initial, 
intermediate-to-long axis ratio deviated noticeably from that of 
the assumed mass model.) We were unable to find any solutions 
with T = 0.75 that did not evolve toward complete axisymme- 
try. We conclude that these solutions do not represent bona-fide 
equilibria. 

Figure ^ shows contours of the projected density of the weak 
cusp models with T = 0.50 and only regular orbits. The con- 
tours remain approximately ellipsoidal until the end of the inte- 
gration. 

5. A MODEL WITH OUTER ANGULAR CONSTRAINTS RELAXED 

We mentioned above that relaxing the angular constraints in 
the outermost shells led to quadratic programming solutions 
with very small errors in the innermost mass cells. However 
these solutions generally exhibited large changes in their shapes 
when evolved forward. We illustrate this in the case of a weak 



cusp nucleus with T = 0.50 and only regular orbits. In the or- 
bital solution for this model, we relaxed the angular constraints 
in the outer 5 shells, i.e. only the masses of these shells were 
fit by the quadratic programming routine, not the individual cell 
masses. By relaxing the outermost angular constraints, we were 
able to find a solution in which log(A) = -2.1 for the innter- 
most shells. Figure || shows the evolution of axis lengths of 
the model. The axis ratios show large fluctuations, especially 
at large radii. Figure ^ shows contours of the projected den- 
sity. At t = 0, the contours are ellipsoidal only at small radii; 
the large-radius contours are irregular due to the large cell mass 
discrepancies in the outer shells. During the integration, this 
model shows significant evolution at large radii, which in turn 
affects the contours at small radii. The final model looks very 
different from the initial model. We found similar behavior in 
other orbital solutions when the outermost constraints were re- 
laxed. 



6. MAXIMIZING THE CONTRIBUTION FROM CHAOTIC ORBITS 

The existence of long-lived triaxial models containing an 
abundance of chaotic orbits is particularly interesting: stars on 
such orbits pass once per crossing time near the center, greatly 
increasing the rate of interactions with the black hole compared 
with spherical or axisymmetric models. We sought to maximize 
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FIG. 7. — Contours of the projected density of the solution with 7 = 1, T = 0.50 and only regular orbits. The upper panels show the contours at t = 0, and the lower 
panels show the contours at the end of the integration. This model was judged stable based on the lack of significant time evolution of the axis ratios (Figure 5); note 
the slight evolution of the contour shapes at the largest radii. 




the chaotic content in our models by minimizing the quantity: 



N M M 

x 2 =^ D i-Y. BiiCi "> 2+ Y. WiCi 



instead of (JT(j). Here W, is the "penalty" associated with the /th 
orbit; increasing Wj tends to decrease the contribution C, from 
the /th orbit in the solution. We set W; = for the chaotic or- 
bits and Wj =W C = (0, 1000, 10000) for the regular orbits. As 
Wc increases, we expect the contribution from chaotic orbits to 





FIG. 9. — Contours of the projected density of the weak cusp solution with T = 0.50 and only regular orbits. The angular constraints of the 5 outer shells are 
ignored. The upper panels show the contours at t = 0, and the lower panels show the contours at the end of the integration. 



increase, although possibly at the expense of the overall qual- 
ity of the fit to the cell masses. We also constructed solutions 
containingonly chaotic orbits. 

Figure and Table [7] show the orbital content for weak and 
strong cusp solutions with T = 0.50 and different values of Wc- 
The bottom panel of Figure [l0| shows solutions containing only 
chaotic orbits. As Wc increases, the regular orbits are replaced 
by high energy (low energy) chaotic orbits in the weak (strong) 
cusp solutions. Evolution of the axis-ratios is shown in Fig- 
ure ]ll]. Remarkably, only the models constructed from purely 
chaotic orbits show substantial evolution in their axis ratios. In- 
spection of the contour plots (e.g. Figure |l2|) does reveal some 
evolution away from elliptical isophotes in the solutions with 
large Wc, but the triaxiality appears robust. We conclude that 
chaotic mass fractions as large as ~ 75% or more might be con- 
sistent with long-lived triaxiality in galactic nuclei. 
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Table 7 

Orbital Content of Solutions with T = 0.5 and 
Different Wc 
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r<0.5 


0.24 


0.06 


0.09 
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0.05 
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0.91 
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0.46 
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0.46 
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0.26 


0.04 


0.04 


0.65 
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0.30 


0.03 


0.03 


0.64 


r < 1.2 


0.30 


0.02 


0.03 


0.65 


7 = 2 : W c = 10000 










r<0.4 


0.15 


0.03 


0.03 


0.79 


r<0.8 


0.18 


0.01 


0.02 


0.78 
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0.17 


0.01 


0.02 


0.79 



scattered into the loss cone, the phase-space region defined by 
orbits with pericenters lying within the black hole's tidal dis- 
ruption radius. In the case of chaotic orbits in a triaxial nu- 
cleus, each passage brings the star near to the center, and the 
time required for a star to pass withi n a distance r, of the black 
hole should scale roughly as r" 1 (e.g. |Gerhard & Binney 1985 ). 
Thus even in the absence of gravitational scattering, the loss 
cone would remain full and the feeding rate could be orders of 
magnitude higher than in axisymmetric nuclei. We will exam- 
ine these "chaotic loss cones" in detail in an upcoming paper 
(Merritt & Poon 2003). 

While our results strengthen the case for triaxiality in galactic 
nuclei, the case for nuclear triaxiality could be made even more 
compelling by the detection of isophotal twists or minor-axis 
rotation at the very centers of galaxies. Such observations will 
be challenging, requiring two-dimensional data on an angular 
scale that resolves the black hole's sphere of influence. Ex- 
isting integral-fi eld spectrographs o n ground-based telescopes 
(e.g. SAURON, |Bacon et al. 2001] ) can only achieve this res- 
olution for the nearest galaxies. Equally valuable would be N- 
body studies demonstrating that triaxial nuclei can form in re- 
alistic mergers. 

M.Y. Poon would like to thank Andrew Mack for stimulating 
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ported by NSF grants AST 96-17088 and AST 00-71099 and 
by NASA grants NAG5-6037 and NAG5-9046. M.Y. Poon is 
grateful to the Croucher Foundation for a postdoctoral fellow- 
ship. 



7. SUMMARY AND DISCUSSION 

We have shown that long-lived triaxial configurations are 
possible for nuclei containing black holes. Models with 
T = 0.5 (maximally triaxial) and T = 0.25 (oblate/triaxial) 
were constructed and found to be stable, retaining their non- 
axisymmetric shapes until the end of the integration interval, 
equal to several crossing times. Models with T = 0.75 (pro- 
late/triaxial) were always found to evolve rapidly to axisymme- 
try; we speculate that prolate/triaxial nuclei do not exist. The 
evolution seen in the nearly-prolate models does not appear to 
be a consequence of orbital chaos; indeed, in our stable solu- 
tions, we were able to replace a surprisingly large fraction of 
the regular orbits by chaotic orbits without inducing noticeable 
evolution in their shapes. We found that at least 50%, and per- 
haps as much as 75%, of the mass could be placed on chaotic or- 
bits in the maximally triaxial and oblate/triaxial solutions. Such 
models violate Jeans's theorem in its standard form (e.g. B in- 
ney & Tremaine 1987) but are consistent with a generalized 
Jeans's theorem ( [Vlerritt 1999| ) if we assume that the chaotic 
building blocks are "fully mixed," that is, that they approxi- 
mate a uniform population of the accessible phase space. This 
appears to be the case for the chaotic orbits in our models which 
have a very short mixing time. While a sudden onset of chaos 
can effectively destroy triaxialit y in models containing a large 
population of regular box orbits (jMerritt & Quinlan 1998j S ell- 
wood 2001), our work shows that at least the central parts of 
galaxies containing black holes can remain triaxial even when 
dominated by chaotic orbits. 

Our results have possibly important implications for the rate 
at which stars are fed to supermassive black holes in galactic 
nuclei. In spherical or axisymmetric nuclei, the feeding rate 
is determined by the rate at which stars on eccentric orbits are 
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FIG. 10. — Cumulative mass fractions F contributed by different kinds of orbits to different shells of the triaxial solutions in which the contribution from chaotic 
orbits has been maximized. The value of Wc is indicated in the upper right hand comer. Left column: T = 0.5, 7 = 1 . Right column: T = 0.5, 7 = 2. 
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FIG . 11. — Evolution of the axis ratios of the orbital solutions illustrated in Figure llOl Columns on the left show models with 7 = 1 , T = 0.50 and various values 
of Wc at r = 1.0 (first column), r = 2.0 (second column), and r = 3.0 (third column). Columns on the right show models with 7 = 2, T = 0.50. 
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FIG. 12. — Contours of the projected density of the weak cusp model with T = 0.50 and Wc = 10000. The upper panels show the contours at t = 0, and the lower 
panels show the contours at the end of the integration. 
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APPENDIX 

EQUIVALENT SPHERICAL MODELS 
We define the equivalent spherical models to have mass density: 

r\-t 



(0=(£) , d 3 = abc. (Al) 



The mass of the central black hole is set to one. The potential is: 

$0) = 27rr--, 7=1 (A2) 
r 

= 4 7 rrf 2 lnf- N )---4 7 r£/ 2 , 7 = 2. (A3) 
\dJ r 

The constant terms in the expressions for the potential were obtained by taking the spherical limits of equation (3) of Paper I. The 
isotropic distribution function f(E) is given by Eddington's formula, 

V2 ( f" d 2 p $ dp . 

1 ' = + hm — (A5) 



4tt 2 \Je d$ 2 \/$-E 

and u = lim,^oo <P(r). We assume that the models extend to infinity. In order to apply Eddington's formula, we need to express p in 
terms of <D. For 7 = 1, we have: 



and 



r ( $ ) = -n ( A6 > 

4-dn 



... d 4d 2 n 
p{§) = - = . (A7) 

r $ + V$ 2 + 8t/7r 



Thus 

d 2 p 4d 2 w 



d<f> 2 ($ 2 + 8c/7r)i 



(A8) 



f(E) = / s r d$. (A9) 

K J E ($ 2 + 8£/7T)i($-£)l 
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Similarly, for 7 = 2: 

1 1 



4ird 2 WO) ' 
1 / $ + 47rc/ 2 (l+lnd/)\ 



r($): 

M= W 6XP V 47Td 2 ) 

and 

d 2 

j0 ($) = — = l6ir 2 d b (W(u)) 2 . 
W(u) is Lambert's W function; it is the inverse of the function u(W) = We w . The distribution function 

dir 2 J E (i+w( M )) 3 V¥^l 

For both models, the differential energy distribution is given by: 



M(E)dE = I6n 2 p{E)f{E)dE, 

/■*"'(£) 

p(E)= / y/2{E-*{r)dr. 
Jo 



